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ABSTRACT 

Results from a direct numerical simulation of laminar flow over 
a flat surface with spherical roughness elements using a spectral- 
element method are given. The numerical simulation approximates 
roughness as a cellular pattern of identical spheres protruding from a 

smooth wall. Periodic boundary conditions on the domain’s horizon- 
tal faces simulate an infinite array of roughness elements extending 
in the streamwise and spanwise directions, implies the parallel-flow 
assumption, and results in a closed domain. A body force, designed to 
yield the horizontal Blasius velocity in the absence of roughness, sus- 
tains the flow. Instabilities above a critical Reynolds number reveal 
negligible oscillations in the recirculation regions behind each sphere 
and in the free stream, high-amplitude oscillations in the layer di- 
rectly above the spheres, and a mean profile with an inflection point 
near the sphere’s crest. The inflection point yields an unstable layer 
above the roughness (where U"(y ) < 0) and a stable region within 
the roughness (where U f, (y ) > 0). Evidently, the instability begins 
when the low-momentum or wake region behind an element, being 
the region most affected by disturbances (purely numerical in this 
case), goes unstable and moves. In incompressible flow with peri- 
odic boundaries, this motion sends disturbances to all regions of the 
domain. In the unstable layer just above the inflection point, the dis- 
turbances grow while being carried downstream with a propagation 
speed equal to the local mean velocity; they do not grow amid the low 
energy region near the roughness patch. The most amplified distur- 
bance eventually arrives at the next roughness element downstream, 
perturbing its wake and inducing a global response at a frequency 
governed by the streamwise spacing between spheres and the mean 
velocity of the most amplified layer. 
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1 Introduction 


Generally, surface roughness promotes transition in the sense 
that under otherwise identical conditions, transition occurs at a lower 
Reynolds number on a rough wall than on a smooth one [Dryden, 1953, 
Sedney, 1973]. Evidently, roughness elements give rise to additional 
disturbances which add to those already present in the laminar, 
boundary layer. A sufficiently rough surface may advance the tran- 
sition location upstream of the smooth-wall position; otherwise, ex- 
tended regions of relatively steady but distorted, laminar flow may 
appear [Sedney, 1973]. Roughness beginning downstream of the loca- 
tion where linear disturbances grow may enhance growth in Tollmien- 
Schlichting ( TS ) waves [Corke & Others, 1986], while disturbances 
from roughness beginning further upstream may “bypass” the known 
linear instability processes [Reshotko & Leventhal, 1981]. 

1.1 Isolated roughness 


The characteristics of flow over three-dimensional roughness ele- 
ments depend on the interaction between the natural boundary-layer 
vorticity and the obstacle [Mason & Morton, 1987]. Upon encoun- 
tering a three-dimensional surface protuberance, the vortex lines may 
concentrate to form a visible vortex core which trails downstream. 
Experiments at low velocity reveal a single smoke streak trailing 
downstream from the element [Mochizuki, 1961]. At higher veloci- 
ties, a horseshoe vortex and a pair of smoke spirals ( chimney vortices 
fixed in space and due to the fluid entering into the rear separation 
pocket) appear close behind the sphere forming two parallel smoke 
filaments trailing downstream parallel to the wall (called the trail- 
ing vortex). As the velocity increases, the trailing vortex filaments 
begin “waving,” while growing over some distance downstream. Fur- 
ther downstream, they cease to grow after oblique stretching by the 
main stream’s shearing stresses. At higher wind speed these vortices 
shed at a dimensionless frequency (non-dimensionalized as a Strouhal 
number, Stk = fk/Uk ) ranging from 0.1 to 0.4 depending on Reyn- 
olds number [Gupta, 1980, Gregory & Walker, 1951]. The trailing 
vortex filaments show waviness at Rek = Ukk/v = 350 and shedding 
at Rek = 700 [Mochizuki, 1961]. (The undisturbed velocity at the 
crest of the roughness equals {/*, and the roughness height equals 
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k . ) In all cases, the dimensionless instability frequency, fu/U^ falls 
above the unstable region of the stability diagram for TS waves. 

l. 2 Distributed roughness 

In contrast to the well understood sequence of instabilities lead- 
ing to transition around isolated roughness, the process for distributed 
roughness remains unknown [Morkovin, 1989]. Although some simi- 
lar mechanisms may exist, the ordered behavior near single elements 
vanishes when neighboring elements are present. We expect an inter- 
action between the streamwise vortices engendered by upstream ele- 
ments and the vortex system of those further downstream; however, 
amalgamation of the wakes from upstream elements should produce 
a low-momentum region amid the roughness patch. This lowers the 
effective roughness Reynolds number, Re £, which depends on the ve- 
locity at the elevation of the roughness in the smooth- wall boundary 
layer. Therefore, with equal roughness Reynolds numbers, a element 
in a distributed roughness patch encounters an extenuated main- 
stream compared to an isolated element. In this case, closely- packed 
elements will not beget the strong horseshoe and chimney vortices 
necessary to initiate the shedding process typically found with iso- 
lated elements until reaching a much higher Reynolds number. 

2 Numerical Simulation of Roughness 

Several numerical studies investigated the mean flow and dis- 
turbance field around isolated three-dimensional roughness elements 
[Mason & Sykes, 1979, Mason & Morton, 1987, Tadjfar, 1990]. Tad- 
jfar employed a three-dimensional triple-deck with time-harmonic 
free-stream disturbances; and Mason et al. solved the full three- 
dimensional, unsteady Navier-Stokes equations. The work presented 
here represents the first computational investigation of distributed 
roughness using both realistic geometry and the three-dimensional, 
unsteady Navier-Stokes equations. 

2.1 Distributed roughness pattern under study 

The distributed roughness under investigation consists of an ar- 
ray of identical spheres protruding from a flat plate according to pre- 
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Figure 1: Side view of the standard density domain showing a 5X5X5 
mesh along the vertical symmetry plane extending from the wall to 
the free stream. 

vious experiments [Kendall, 1981,Tadjfar & Others, 1985]. Whereas 
the experimental studies consisted of a spatial array of a large num- 
ber of spheres, the computational domain consists of the single cell 
shown in figures 1 and 2. 

2.2 Boundary conditions 


All previous experiments with distributed roughness were done 
in spatially-developing boundary layers or “open” flows. By applying 
periodic boundary conditions to the horizontal faces of a single cell 
and an impermeable boundary condition along the cell’s top, we 
simulate an infinite array of roughness elements extending upstream, 
downstream, and spanwise. These boundary conditions effectively 
“close” the domain and result in a parallel flow. 

In simple geometries, periodic conditions allow us to simulate 
spatial development by interpreting the time-dependent data as that 
coming from a domain moving downstream with the group velocity 
of the disturbances. This “moving box” approach is often effective 
in these geometries [Singer & Others, 1986, Spalart & Yang, 1987, 
Laurien & Kleiser, 1989]. 

This artifice fails in complex geometries, however. The domain 
must remain fixed at a given streamwise position, and temporal data 



Figure 2: Plan view of the standard density domain showing a 
5X5X5 mesh along a horizontal plane passing through the center 
of the spheres. 

can only be used to indicate frequency and phase content along with 
relative amplitude between signals from different locations; it cannot 
reveal anything about growth rates since disturbances cannot grow at 
a fixed location in a flow with steady mean quantities. The state with 
steady mean flow is called the time-asymptotic state. This restriction 
does, however, provide a convenient check on the numerical scheme: 
After achieving the time-asymptotic stage, a solution showing grow- 
ing fluctuations indicates poor spatial or temporal resolution, while 
decaying fluctuations cannot occur. (These comments do not apply 
to time-asymptotic flows with imposed disturbances. In these cases, 
a time-asymptotic flow may be “perturbed” by introducing artificial 
disturbances which may grow or decay in time.) 

Finally, a word about roughness and TS waves is appropriate. 
Although experiments reveal enhanced growth in the TS band of 
frequencies when roughness begins downstream of the critical point 
for TS wave growth, the periodic and free-stream boundary condi- 
tions applied to the single roughness cell shown in figure 1 preclude 
TS waves — TS waves have wavelengths which typically equal ten 
boundary-layer thicknesses and amplitudes which decay asymptoti- 
cally with distance from the wall. Unfortunately, the finite height do- 
main, imposed free-stream boundary conditions, and relatively short 
cell length are incongruous with this type wave. 
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2o3 Numerical procedure 


Both the incompressible continuity equation, given by 



and the incompressible Navier-Stokes equation, given by 

du{ du{ _ 1 dp d 2 Ui 

dt + dxj p dx{ V dxjdxj ’ 


( 1 ) 

( 2 ) 


are solved in the domain described in sections 2.1 and 2.2 using 
a spectral-element code. The spectral-element method employs a 
three-step, time-splitting scheme where the non-linear, pressure, and 
viscous terms in the incompressible Navier-Stokes equations are writ- 
ten in separate fractional steps [Orszag & Kells, 1980]. By introduc- 
ing intermediate velocities, these steps yield consecutive solutions 
based first on the non-linear, then on the pressure, and finally on the 
viscous terms. 


2.4 Input parameters 


When the incompressible Navier-Stokes equations (see equa- 
tion 2) are written in non-dimensional form, three dimensionless 
quantities appear: the Reynolds number, Re = UL/ir, the Strouhal 
number, St = fL/U ; and the Froude number, F r = U/y/gZ . With 
periodic boundary conditions, a body force generates the velocity, 
and a relationship exists between the Reynolds and Froude num- 
bers. Furthermore, in a flow without external forcing, the Strouhal 
number governs unsteadyness due to naturally arising instabilities. 
Hence, the Reynolds number or, for a given geometry, the fluid vis- 
cosity and the body force, completely determines the flow. 

Here we employ a special body force which, when twice differ- 
entiated with respect to y, yields the Blasius profile in the absence 
of roughness elements. It is derived by considering the streamwise 
momentum equation in a boundary layer on a smooth, flat surface 
with no imposed streamwise pressure gradient. Equation 2 in this 
case reduces to 


du du d 2 u 

U te +V d-y =U W 


( 3 ) 
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By interpreting the right-hand side of equation 3 as a body force, we 


get 


du du 

u ~+ v di = /*(»)» 


dx 


( 4 ) 


with 


f B (y) = v 


d^UBjy) 

d 2 y 


( 5 ) 


while the subscript B refers to the Blasius value. By twice differenti- 
ating the Blasius velocity profile according to equation 5, we obtain 
a body force which yields the Blasius profile on a smooth, flat surface 
with no imposed pressure gradient. When this body force is applied 
to a domain with surface roughness elements, the resulting velocity 
field represents the influence of the roughness elements on the Blasius 
boundary layer. 


2.5 Time asymptotic state 


From a given initial state, a certain time elapses while mean- 
flow quantities develop towards a time-asymptotic state. The time 
required to reach this state depends on the starting solution as well 
as the longer of the convection and diffusion times. 

The convection time to reach the asymptotic state depends on the 
time for a particle to pass across the domain, given by r c = L x /Uoo, 
where Uqo equals the streamwise free-stream velocity, and L x equals 
a streamwise length scale. 

The diffusion time to reach the time-asymptotic state depends 
on momentum diffusion across the layer. For a fixed boundary-layer 
thickness and viscosity, momentum diffusion across a boundary layer 
of thickness 6 requires 6 2 fv time. 

The ratio of diffusion to convection times equals 


Td _ <S 2 /^ 

T c WC/oo 


<£> , 


U„S / - Xp 

= (-r)Res- 


_6_ 

L x 


( 6 ) 


Typically, the Reynolds number equals several hundred, and the 
ratio of boundary-layer thickness to streamwise length is of order 
one; hence, diffusion requires two orders-of-magnitude more time 
than convection. For typical flows of interest, diffusion requires 
10 3 time units. Since the first-order time-accurate spectral-element 
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Figure 3: Side view of the high-density domain showing a 5X5X5 
mesh along the vertical symmetry plane extending from the wall to 
the free stream. 

method requires small time steps — typically At < 0.01 — for good ac- 
curacy, approximately 10 5 time steps elapse before reaching a time- 
asymptotic flow. 

2.6 Cases studied 


Some previous experiments using distributed roughness elements 
[Gartshore & de Croos, 1979, Lee & Soliman, 1977, Kendall, 1981] 
found that roughness “density,” or spacing between elements, con- 
trols certain mechanisms responsible for transition. Based on this 
idea, I selected three different roughness densities for study: stan- 
dard, low, and high. The roughness density was varied by stretching 
or compressing the domain in the streamwise direction while holding 
the spanwise dimension, the number of finite elements, and the to- 
tal number of grid points constant. Consequently, the high-density 
domain was shorter (see figure 3), and the low-density domain was 
longer (see figure 4) than the standard domain. The standard-density 
domain consists of a square with streamwise and spanwise lengths of 
8.4; the high-density domain consists of a rectangle with a streamwise 
length of 6.0 and spanwise length 8.4; and the low-density domain 
has a streamwise length of 12.0 and a spanwise length of 8.4. All 
cases used Re k = 175, k/6* = 0.72, and the Blasius based body force 
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Figure 4: Side view of the low-density domain showing a 5X5X5 
mesh along the vertical symmetry plane extending from the wall to 
the free stream. 

(see section 2.4) to drive the flow. For each case, the one-cell grid 
consisted of 144 elements and a 5X5X5 mesh within each element. 
The standard-density case was also run with a finer, 7X7X7, mesh. 

The calculations performed on a single-cell domain correspond to 
the assumption of “strong periodicity” [Karniadakis & Others, 1988]. 
Unfortunately, we have no a priori knowledge whether or not a given 
flow will behave differently when the number of cells are changed. To 
gain insight into the effect of the streamwise, periodic boundary con- 
ditions on the nature of the disturbance wavelength and frequency, a 
two-cell high-density case and a two-cell low-density case were run. 
These two-cell cases contained two 144 -element grids aligned in the 
streamwise direction, used the 5X5X5 mesh, and were twice as long 
as their single-cell counterparts. A change in results between the 
one- and two-cell cases will indicate the presence of a sub-harmonic 
disturbance. 

3 Steady State Results 


Since we compute the unsteady Navier-Stokes equations until 
a time-asymptotic state is reached, all of the results are unsteady 
above a critical Reynolds number. To distinguish steady results from 
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unsteady results (given in the next section) the term u steady state” 
is used loosely to refer to quantities which do not vary significantly 
with time. These include mean flow boundary-layer quantities and 
instantaneous velocity and vorticity contours. 

3.1 Mean flow results 


Results for the three different roughness densities using a 5X5X5 
mesh with one cell are given in table 1. Also appearing in table 1 
are results for the standard-density case using a 7X7X7 mesh. The 
corresponding smooth- wall, Blasius results are given for comparison. 
Some of the quantities appearing in the table need to be defined be- 
fore discusing the results. The ratio of rough-wall to smooth-wall 
friction coefficients equals . The maximum velocity divergence, 
equals ( V- V) m , and the velocity divergence integrated over the entire 
domain, equals (V • V ) a — both indicate the quality of the numerical 
solution, since they equal zero in incompressible flow. Finally, the 
displacement and momentum thicknesses appear here in dimension- 
less form. 

According to table 1, an improvement in grid resolution from 
a 5X5X5 to a 7X7X7 mesh per element yields a slight increase in 
drag, evidence that the course mesh fails to capture the smaller scales 
required to resolve the flow. The results also show an increase in 6* 
and a decrease in 6 compared to the smooth-wall (Blasius) values 
and consistent with inflected or adverse pressure-gradient profiles. 


density 

mesh 

S m 

e 

H 

(V • V) m 

(V-V) a 

~TW~ 

(Cth 

stand. 

5X5X5 

1.85 

.508 

3.64 

0.0200 

-0.0176 

2.62 

stand. 

7X7X7 

1.88 

.530 

3.57 

— 

— 

2.64 

high 

5X5X5 

1.80 

.473 

3.80 

0.0148 

-0.0128 

2.69 

low 

5X5X5 

1.84 

.529 

3.48 

0.0322 

-0.0352 

2.49 

Bias. 

theory 

1.72 

.664 

2.59 

0.0 

0.0 

1.0 


Table 1: Mean flow results for the various roughness densities 
and the corresponding smooth-wall, Blasius data 


Compared to the standard domain, the high-density case effects 
an increase in ^y 3 and H and a decrease in both 6* and 0, the 
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Figure 5: Side view of the streamwise velocity contours on a 

standard-density domain with a 5X5X5 mesh along the vertical sym- 
metry plane extending from the wall to the free stream. 

opposite occurs in the low-density case. 


3.2 Velocity contours standard density domain 


The standard-density domain velocity contours are given in fig- 
ures 5 through 10; vorticity contours are given in figure 11. Figures 5 
and 6 show the streamwise and vertical velocity on a plane through 
the center of the middle sphere. A small reverse flow appears in the 
region behind the sphere. The vertical velocity (figure 6) equals zero 
above two roughness heights while remaining within 4.4 percent of 
the free-stream velocity near the elements. The small upward veloc- 
ity behind the element remains confined to the region between the 
centerline and crest. At this Reynolds number, the sign and location 
accede with the chimney vortex of an isolated roughness element — 
a vortex rotating counter to and above the horseshoe vortex. In 
this case, however, the neighboring roughness elements beget a low- 
momentum region and prevent the formation of robust horseshoe and 
chimney vortices. Instead, a flow dominated by kinematics results: 
continuity in concert with the solid-wall boundary conditions. The 
upward velocity behind the sphere depends on continuity: the fluid 
passing around the sides of the element converges in the wake and 
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Figure 6: Side view of the vertical velocity contours on a standard- 
density domain with a 5X5X5 mesh along the vertical symmetry 
plane extending from the wall to the free stream. 

gets forced upward (the wall prevents it from going down). 

Reverse-flow regions indicate the importance of inertia since no 
reverse-flow regions exist in a viscous dominated flow. Thus, the 
upward and upstream velocity in the wake indicates that a chimney 
vortex, albeit a weak one, forms. Does the horseshoe vortex ap- 
pear? With an isolated sphere, the horseshoe vortex begins in the 
stagnation region and extends into the wake. In distributed rough- 
ness, however, other roughness elements disturb the velocity in the 
wake; consequently, the stagnation region becomes the best place to 
evince the horseshoe vortex. Unfortunately, only an extremely weak 
horseshoe vortex exists. A true horseshoe vortex would appear in the 
stagnation region as a downward velocity near the element’s center- 
line, an upstream velocity near the wall, and an upward velocity out 
in front of the element. As shown in figure 6, the downward velocity 
occurs near the sphere; however, neither the upward nor upstream 
velocities are evident. 

The same velocity components on a horizontal plane passing 
through the center of the spheres are shown in figures 7 and 8. A high 
streamwise velocity ‘‘channel” appears between the rows of spheres 
as a consequence of forcing the flow along the geometric symmetry 
plane, and a reverse-flow region appears downstream of each ele- 
ment. The vertical velocity contours show downward flow ahead and 


12 



-0,0025 
0J0025 
OjOC75 
OX) 125 
0X3175 
0X3225 
0X3275 
0X3925 
01)975 
0XM25 
ODC75 
OXX52S 
0X3875 
OD625 
0XX75 
0X3725 
0X3775 
0X3825 
0X3875 


Figure 7: Plan view of the streamwise velocity contours on a 

standard-density domain with a 5X5X5 mesh along a horizontal 
plane passing through the center of the spheres. 
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Figure 8: Plan view of the vertical velocity contours on a standard- 
density domain with a 5X5X5 mesh along a horizontal plane passing 
through the center of the spheres. 
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Figure 9: Cross-sectional view of the vertical velocity contours on a 
standard-density domain with a 5X5X5 mesh along a vertical plane 
passing through the domain’s inlet extending from the wall to the 
free stream. 

upward flow behind each element, consistent with the continuity ar- 
gument previously outlined. 

Figures 9 and 10 present the vertical and spanwise velocity on a 
cross-sectional plane along the inlet of the domain. This plane cuts 
through the center of the corner spheres and lies just downstream 
of the center sphere in an imaginary cell upstream of the computa- 
tional cell. Hence, this plane shows the central sphere’s wake struc- 
ture. The spanwise velocity always points toward the centerline, 
clearly showing how the fluid fills in the region behind the center 
sphere while being pushed inward by the corner spheres — both ef- 
fects act in the same direction, yielding a strong spanwise velocity. 
The vertical velocity pattern illicit s more interest. The central re- 
gion shows upward velocity consistent with the continuity constraint 
discussed previously. A larger upward velocity region develops above 
each sphere — a manifestation of the inertia of the fluid passing over 
the obstacle. (If the flow were free of inertia, all velocity on this plane 
would be in the streamwise direction.) The two downward velocity 
regions reflect wake filling from above the corner spheres. 

Streamwise vorticity on a cross-sectional plane at the inlet is 
shown in figure 11. These vorticity regions seem to match those pro- 
duced by the horseshoe and chimney vortices behind isolated rough- 
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Figure 10: Cross-sectional view of the spanwise velocity contours 
on a standard-density domain with a 5X5X5 mesh along a vertical 
plane passing through the domain’s inlet extending from the wall to 
the free stream. 
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Figure 11: Cross-sectional view of the streamwise vorticity contours 
on a standard-density domain with a 5X5X5 mesh along a vertical 
plane passing through the domain’s inlet extending from the wall to 
the free stream. 
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Figure 12: Side view of the vertical velocity contours on a high- 
density domain with a 5X5X5 mesh along the vertical symmetry 
plane extending from the wall to the free stream. 

ness elements. Congruent with the vorticity system downstream of 
the center sphere without neighbors, on the left we see positive vor- 
ticity below and negative vorticity above the midpoint; the signs 
reverse on the right. We could easily link the trailing vortices from 
the sphere upstream and the vorticity pattern on the inlet plane. 
However, if these regions actually contain vortices, circular velocities 
should appear in the spanwise and vertical contour plots. We see no 
evidence of circular motion. 

3.3 Velocity contours high and low density domains 

The vertical velocity contours on the symmetry plane (see fig- 
ures 12 and 13 ) show that the dense-packing maximum velocity 
decreases and the sparse- packing value increases compared to the 
standard domain. Evidently the sparse array allows an increase in 
momentum between the wake of one sphere and the stagnation region 
of the next. Further confirmation appears in the streamwise velocity 
contours on the horizontal plane passing through the sphere cen- 
ters (not shown). The maximum streamwise velocity equals 0.0875 
in the standard domain, 0.054 in the dense domain, and 0.115 in 
the sparse domain, a difference of a factor of two between the dense 
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Figure 13: Side view of the vertical velocity contours on a low-density 
domain with a 5X5X5 mesh along the vertical symmetry plane ex- 
tending from the wall to the free stream. 

and sparse cases. Clearly the dense array creates a large region of 
low-momentum fluid amid the spheres. An even more striking ef- 
fect is evident in the spanwise and vertical velocity contours on the 
inlet cross-sectional plane (figures 14 and 15). Figure 14 shows the 
spanwise velocity contours for the dense array. The usual conver- 
gent velocity regions are shown; however, an outward velocity below 
this region exists, a downward velocity between these regions ap- 
pears, and an upward velocity outside the downward velocity region 
develops — the horseshoe vortex. There is no evidence of the chim- 
ney vortex. Somehow the close proximity of the corner spheres to 
the wake of the center sphere generated the horseshoe vortex. 


4 Time Dependent Results 


During each run, we store streamwise, spanwise, and vertical 
velocity signals from five numerical probes. Probe locations are given 
in figures 16 and 17. Probe number 48, located behind the center 
sphere, lies within one sphere elevation from the wall; probe numbers 
62, 78, and 88 lie approximately one and one-half sphere diameters 
from the wall; and probe number 120 lies two and one-half diameters 
above the wall. Probe numbers 62 and 78 occupy symmetric positions 
behind each of the two forward corner spheres. Similarly, probes 78 
and 88 occupy equivalent positions behind the center sphere and one 
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Figure 14: Cross-sectional view of the spanwise velocity contours on 
a high-density domain with a 5X5X5 mesh along a vertical plane 
passing through the domain’s inlet extending from the wall to the 
free stream. 
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Figure 15: Cross-sectional view of the vertical velocity contours on 
a high-density domain with a 5X5X5 mesh along a vertical plane 
passing through the domain’s inlet extending from the wall to the 
free stream. 
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Figure 16: Side view of the probe locations 



Figure 17: Plan view of the probe locations 
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Figure 18: Vertical (top figure), streamwise (middle figure), and 
spanwise (bottom figure) velocity versus time for the five probes in 
the standard-density domain. 

of the forward corner spheres. Thus, differences in the signals from 
each of these probes may reveal flow asymmetries and phase shifts 
in the disturbance waves. 

4.1 Standard density domain 

Using the standard domain, mesh dependent solutions revealed 
both steady and unsteady time-asymptotic states. Unsteady flows 
using a 5X5X5 mesh decayed with a 7X7X7 mesh. 

Results using a 5X5X5 mesh are given in figure 18. Although 
not shown, the total kinetic energy variation (an indication of flow 
development toward the time- asymptotic state) was less than 0.1 
percent over the time contained in the graphs. This corresponds to 
a dimensionless convection time, r c = AtUoo/ L x , of 66, or a dimen- 
sionless diffusion time, = A tv/k 2 , equal to 0.51. These numbers 
indicate that a fluid element in the free stream passed across the 
domain 66 times, but momentum only diffused across approximately 
0.7A: during the measuring interval. The small change in total kinetic 
energy, however, indicates that the flow may appropriately classified 
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Figure 19: Vertical velocity amplitude for each of the five probes 
versus Strouhal number in the standard- (top figure), high- (middle 
figure), and low-density (bottom figure) domains. 

as time-asymptotic. No oscillations developed in any of the velocity 
components from probe numbers 48 and 120 until the r.m.s. value of 
probe 78’s streamwise velocity reached 0.5 percent of the free-stream 
velocity. At this point probe 120 began a low-level oscillation. Ver- 
tical, streamwise, and spanwise oscillations with a 180° phase shift 
occurred in probes 78 and 88; a small phase shift exists between 
probes 62 and 88. These oscillations produce a strong single spike at 
Stk = kf /Uk = 0.33 in the standard domain’s frequency spectrum 
(see figure 19). This corresponds to the “passing” frequency, Uk/ L x , 
for the time a parcel of fluid at elevation y = k requires to cross 
between spheres spaced at a distance L x with a Blasius velocity Uk . 

4.2 High and low density domains 

If the roughness density controls the characteristic frequency, a 
lower oscillation frequency would exist in a low-density domain, and 
a higher one in a high-density domain. The high-density domain’s 
streamwise length equals 6.0; the low-density domain’s streamwise 
length equals 12.0. These lengths correspond to passing frequency 
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Figure 20: Vertical (top figure), streamwise (middle figure), and 
spanwise (bottom figure) velocity versus time for the five probes in 
the high-density domain. 

based Strouhal numbers of 0.47 and 0.23 for the high- and low-density 
domains respectively. The time-dependent data from the five probes 
are given in figures 20 and 21 for these cases. Although not obvious 
from these figures, the dominant oscillation frequency, according to 
the frequency spectra shown in figure 19, does agree with the pre- 
dicted Strouhal numbers given above. 

5 Analysis 

5.1 Mean flow variations with roughness density 

At first glance it would appear that the displacement and mo- 
mentum thicknesses would increase as the roughness elements be- 
come more densely packed. However, in the limit of infinite roughness 
density, these quantities would equal their Blasius values since the 
roughness effectively becomes a flat plate. A plot of 6 * and 6 versus 
roughness density would show Blasius values at both the infinite- 
density and zero-density limits. In between, the displacement thick- 
ness increases and the momentum thickness decreases compared to 
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Figure 21: Vertical (top figure), streamwise (middle figure), and 
spanwise (bottom figure) velocity versus time for the five probes in 
the low-density domain. 


Blasius values, consistent with inflected or adverse pressure gradient 
profiles. 

Another parameter which changes with roughness density is the 
friction coefficient. Schlichting investigated how the flow resistance 
depends on roughness density [Schlichting, 1936]. He defined rough- 
ness density as 


F r = 


projected area of roughness on plate 
plate area 


( 7 ) 


and found that the maximum resistance does not occur at the max- 
imum roughness density but at a considerably lower value: 


(Fr)max — 0.4. 


(8) 


The three roughness densities used in this study correspond to F r = 
0.20 in the standard domain, F r = 0.28 in the high-density domain, 
and F r = 0.14 in the low-density domain. Schlichting defines a re- 
sistance coefficient in terms of the additional drag due solely to the 
roughness elements as W T and the corresponding friction coefficient 


due to the roughness elements as 


W r 


(Cf)roughonly — 1^2^ ’ 


( 9 ) 


Defining the ratio of rough- to smooth- wall friction factors as \ gives 

( r* \ (Cf)smooth(X “ 1) / in \ 

\^f /rough only = J • (10) 

Using the given values of F r and the results from table 1 for x gives 
the results shown in table 2. Also included are Schlichting’s experi- 
mental results. 


Case 

F r 

Uk/U CO 

(C/ trough only 

standard-density 

0.20 

0.43 

0.019 

high-density 

0.28 

0.43 

0.013 

low-density 

0.14 

0.43 

0.023 

Schlichting 

0.126 

0.704 

0.014 

Schlichting 

0.349 

0.678 

0.017 

Schlichting 

0.907 

0.826 

0.007 

Schlichting 

0.126 

0.570 

0.011 


Table 2: Roughness friction coefficient versus roughness density, 
comparison between numerical results and experiments. 


Although the smallest total friction factor occurs in the low- 
density domain, we see that the friction factor due to the spheres 
alone is greatest in this case. The low-density domain permits more 
acceleration between spheres resulting in a faster approaching stream; 
consequently, each element produces greater resistance. 

The numerical cases studied laminar flow and yielded rough- 
ness resistance coefficients higher than Schlichting’s turbulent-flow 
results. This comports with expected behavior: an increase in fric- 
tion coefficient with decreasing Reynolds number. 

5.2 Rough to smooth wall comparison 

The previous discussion concerns the effect of roughness den- 
sity on the rough-wall boundary-layer integral quantities; this sec- 
tion relates the roughness quantities to the corresponding Blasius or 
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smooth- wall values. All cases reveal a larger dimensionless displace- 
ment thickness than the Blasius value 1.72, a dimensionless momen- 
tum thickness smaller than the Blasius value 0.664, and a shape 
factor larger than the Blasius value 2.59. These results reflect the 
nature of the roughness boundary layer: a low-momentum region 
surrounding the roughness elements and a high-momentum region 
just above the roughness. Evidently, the spheres push fluid toward 
the free stream, leaving the region within one diameter from the wall 
with low momentum. Near the wall, the term (1 — U/Uoo ) is large, 
yielding a larger 6*. Also near the wall, U/Uoo is small, yielding a 
lower 0 . 

The difference between the rough- and smooth- wall displacement 
thicknesses indicates the extent of the velocity profile displacement 
due to roughness. For the smooth wall we can estimate the dis- 
placement thickness as 6* « £/3 = 3.67; a typical rough- wall mea- 
sured value gives 6* « 3.90. Therefore, the additional displace- 
ment due to roughness normalized by the roughness height equals 
A S*/k = (3.90-3.67)/2.8 = 0.082. Comparing this measured change 
in displacement thickness to the volumetric average thickness of the 
roughness itself, 


Ah. _ _ 

k kL x L z 


( 11 ) 


shows that, roughness displaces the Blasius layer an additional three 
quarters of the average roughness height, AS* /Ah « 0.75. 


5.3 Streamwise vorticity 

The experiments with isolated roughness show the importance of 
trailing vortices on boundary-layer instabilities. With this in mind, it 
is of interest to examine the distributed roughness flow for streamwise 
vortices. With the aid of figures 9 and 10, the spanwise and verti- 
cal velocity contours show no evidence of circular flow. Instead, the 
vorticity patterns depend on the gradients of mean velocity. Mathe- 
matically, the streamwise vorticity component equals 

C x = -( dw/dy - dv/dz). (12) 

In the results with spherical roughness, the variation of spanwise 
velocity with vertical elevation, dw/dy , dominates the other term. 
Referring to figure 10, we see that w grows in magnitude from zero 
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at the wall to a maximum at the sphere’s midplane and back to zero 
above two diameters. On the right side of the plane, w is positive; 
hence, dw/dy > 0 near the wall, and the vorticity is negative. Above 
the midplane, dw/dy < 0, and the vorticity is positive. The opposite 
signs correspond to the left half of the plane. This pattern exactly 
matches the streamwise vorticity regions of figure 11. Evidently, the 
presence of other spheres and the relatively low Reynolds number 
used in the computations prevent the strong streamwise vortices ob- 
served in experiments with isolated roughness. 

5.4 Oscillation frequency and phase 


The instability scenario begins when the local recirculation re- 
gions behind each sphere become unsteady above a certain critical 
Reynolds number (still unknown). This unsteady motion affects the 
mean flow by introducing pressure and velocity disturbances which, 
due to the periodic boundary conditions and incompressibility of the 
fluid, instantly fill the entire domain. 

Once introduced, these disturbances grow or decay while prop- 
agating downstream with the local convection velocity of the mean 
flow. Eventually, those disturbances receiving the most energy in 
unstable modes dominate the flow. It appears that growth does not 
occur in the recirculation regions behind each sphere, since there 
is little incentive for amplification this close to the wall. Instead, 
the major growth occurs in the destabilized region just above the 
inflection point. 

This comports with the stability characteristics governed by U"(y). 
Inflected profiles usually have a higher amplification rate for the same 
frequency and Reynolds number than a profile with purely negative 
U"(y) [Obremski & Others, 1969]. Although inflected profiles are 
less stable, there is a stabilizing effect in the region between the in- 
flection point and wall where £/"(y) >0. An examination of the 
streamwise momentum equation reveals that U"(y) acts like a force 
per unit mass on a parcel of fluid. When positive, as within the 
region between the wall and the inflection point, it accelerates or 
stabilizes a fluid parcel; outside of the inflection point it opposes 
or destabilizes the motion. This explains why oscillations at probe 
48, located below the inflection point, are small; why oscillations at 
probes 62, 78, and 88, located just above the inflection point, are 
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large; and why oscillations at probe 120, located in a region where 
U"(y) is small, are small. What explains the oscillation frequency? 

Apparently, the most amplified disturbances occur in the layer 
just above the spheres, presumably carried downstream with a ve- 
locity of the local mean flow. Upon reaching the exit boundary, the 
disturbances reappear along the inlet boundary then resume their 
journey downstream. Once again, they encounter the exit bound- 
ary, and the pattern repeats ad infinitum. A disturbance interacts 
with the low-momentum recirculation zones behind a sphere each 
time it passes overhead. (It also interacts with higher momentum 
regions, but these regions do not respond significantly.) If the flow 
is unstable, the interaction increases the amplitude of the distur- 
bance and causes the recirculation zone to move more vigorously. 
This generates a periodic motion with period equal to the time for 
a fluid element to pass between spheres at the elevation of the most 
amplified disturbance — in this case the period equals L x /Uk- An in- 
dividual recirculation region interprets this as a periodic forcing or 
“buffeting.” If this buffeting frequency falls near the body’s natural 
frequency, a response at the buffeting frequency occurs. This is re- 
ferred to as the “lock-in” mode [Karniadakis & Triantafyllou, 1989]. 
When the buffeting frequency differs greatly from the natural fre- 
quency, the flow responds at the isolated-body natural frequency, 
referred to as “non-lock-in.” Evidently, the three roughness density 
cases under consideration — standard, high, and low — all satisfy the 
lock-in mode criteria. In each case, the oscillation frequency equals 
the forcing frequency. 

Presumably, a higher or lower density case would show the non- 
lock-in mode. In this case, the oscillation frequency would equal 
the isolated element natural frequency and would be independent 
of periodic boundary conditions or element spacing. We may also 
expect a shear-layer based instability instead of one based on vortex 
shedding. 

According to Karniadikis and Triantafyllou, 1989, the distur- 
bance or forcing amplitude also influences the response. They present 
a plot of forcing amplitude versus forcing frequency, with unsteady 
regions labeled “lock-in” and “non-lock-in.” The lock-in region oc- 
curs within an open-end-up parabolic curve centered at the natural 
frequency of the isolated body. The non-lock-in region surrounds the 
lock-in zone. Furthermore, the unsteady region — whether lock-in or 
non-lock-in — only occurs above a threshold amplitude. Thus, forcing 


27 


at the natural frequency with an amplitude just above the threshold 
will produce a lock-in response; forcing with the same amplitude at 
a slightly different frequency will produce a non-lock-in response. As 
the forcing amplitude increases, the lock-in mode region grows; pre- 
sumably, at infinite amplitude any forcing frequency will yield the 
lock-in behavior. 

In the numerical calculations shown here, the only forcing is due 
to numerical approximations, and it seems reasonable to assume that 
the forcing level is greater with a 5X5X5 mesh than with a 7X7X7 
mesh. Therefore, the lock-in region is larger with the 5X5X5 mesh 
(possibly the lack of any response with the 7X7X7 mesh indicates 
that the disturbance level with a 7X7X7 mesh falls below the thresh- 
old amplitude). We conclude that the unsteady results will always 
depend on the forcing amplitude or mesh resolution in flows with no 
externally applied disturbances. 

We’ve discused the oscillation amplitude and frequency, can we 
learn anything from the phase between different probes? In the high- 
density domain (see figure 20), probes 62, 78, and 88 oscillate in all 
three directions while probes 48 and 120 show weak or no oscillations. 
In all directions, probes 62 and 78 oscillate with approximately a 
180° phase shift, and probes 78 and 88 oscillate with a 90° phase 
shift. The 90° phase shift occurs at a frequency of 0.054 cycles per 
time corresponding to a time shift of 4.6 units between peaks. The 
streamwise distance between probes 78 and 88 equals approximately 
3.0 length units in the high-density domain; therefore, the streamwise 
phase speed equals 

c * = IS! = °- 65 - (i3) 

The Blasius velocity at the elevation of these probes equals 0.62, 
further evidence that the disturbances propagate with the local con- 
vection velocity. A disturbance with this frequency and phase speed 
corresponds to a wavelength of 

A = ^ = 12.1. (14) 

This wavelength equals twice the domain length, and the 90° phase 
shift between probes 3.0 length units apart corresponds to one quar- 
ter of a wave with a wavelength of 12.0 length units. 

How can a disturbance with wavelength equal to twice the do- 
main length exist in a domain with periodic boundary conditions? 
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As in all numerical simulations, disturbances must satisfy the bound- 
ary conditions. This may be a restriction in one-dimensional com- 
putations; however, periodic boundaries may not restrict the dis- 
turbance character in two- or three-dimensional domains when the 
disturbances display only narrow spanwise or vertical extent. When 
a small disturbance travels downstream at an angle oblique with re- 
spect to the direction of periodic boundary conditions, or when it 
rises toward the free stream while moving downstream, there’s no 
need for the disturbance to match inlet boundary conditions upon 
reaching the exit — it moved relative to its original position. Hence, 
the disturbance wavelength may be longer or shorter than the do- 
main’s length. Evidently, this is what occurs with small roughness 
elements. 

The oscillations change in the low-density domain: probes 48 and 
120 oscillate more vigorously than in the high-density domain, and 
the phase shifts between different probes no longer remain equal in 
the three directions (see figure 21). In the vertical direction, the 
phase shift between probes 88 and 78 approximately equals 5° (cor- 
responding to 4 time units); the phase shift between probes 62 and 
78 equals 10° (corresponding to 8 time units). Although probe 48’s 
results are not shown, probes 48 and 120 oscillate nearly in-phase 
with a 5° to 10° phase shift from the other three. A zero phase shift 
exists between probes 62, 78, and 88 in the streamwise direction. In 
the spanwise direction, a small 5° phase shift exists between probes 
78 and 88 and nearly a 180° phase shift between these probes and 
probe 62. Clearly the simple structure evident in the high-density 
case is lost when the spheres are spaced further apart. 

Can we classify the various waves appearing in the different cases? 
Instabilities in periodic domains involve either standing or traveling 
waves. Since the standing wave’s phase speed equals zero (allowing 
independence of frequency and wavelength) a flow field dominated by 
a single standing wave reveals signals from probes located at different 
points with either a 0° or 180° phase shift. On the other hand, 
a traveling wave must must satisfy = / A, where the frequency 
equals / and the wavelength equals A. In this case, any phase shift 
between probes is possible. 

With the periodic pattern of spheres, the recirculation zone be- 
hind each element at a Reynolds number below that associated with 
vortex shedding contains a standing wave: the entire separation zone 
oscillates at a single frequency. This motion influences those particles 
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outside the recirculation zone, but, since they fall in a region domi- 
nated by convection, the standing- wave description would not apply, 
and the flow contains a traveling wave. A probe, however, measures 
the same frequency within the recirculation zone or outside in the 
convection zone; although, the two signals differ in phase — phase 
speed equals zero within the recirculation zone, and is finite outside. 

The vortex-shedding mode at higher velocity contains no stand- 
ing waves. In this regime, shed vortices convect downstream, and 
numerical probes located at different positions in the wake measure 
a phase shift between signals corresponding to the time for shed 
structures to pass. 

In the high-density domain the 180° phase shift between probes 
62 and 78 indicate a standing wave in the horizontal direction, while 
the 90° shift between probes 78 and 88 indicates a traveling wave in 
the streamwise direction. 

In the low-density domain nearly all the signals from the vari- 
ous probes show small phase shifts corresponding to traveling waves. 
Only the spanwise signal from probes 62 and 88 show the character- 
istic 180° phase shift associated with standing waves. The non-zero 
but small phase shifts correspond to a rapid traveling wave; the 180° 
shift may indicate a standing wave or a traveling wave with wave- 
length equal to the domain length. 

Apparently, the high-density domain’s disturbances behave like 
standing waves, and the low-density domain’s disturbances behave 
more like traveling waves. This indicates a prominence of vortex 
shedding in the low-density domain and vortex waving in the high- 
density domain. 

5.5 Two cell behavior 


What happens to the characteristic frequency after adding an- 
other cell to the domain? Two possibilities exist: the oscillation 
frequency drops by a factor of two, as a sub-harmonic oscillation 
with wavelength equal to twice the one-cell length occurs; or the 
frequency remains the same. 

In the first scenario, the periodic boundary conditions must in- 
fluence the flow field, changing the number of cells probably changes 
the solution, and the individual roughness elements do not directly 
produce disturbances leading to oscillations; instead, a shear layer 
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based disturbance takes place. In the second scenario, it appears 
that each roughness element contributes to the disturbance, the os- 
cillation seems directly linked to the spacing between elements, and 
no solution variations occur with changes in cell number. 

Although the two-cell results are not shown, the two-cell high- 
density domain yields St = 0.49, and the two-cell low-density do- 
main yields Stk = 0.24. These agree with the one-cell results, and 
it appears that the single cell adequately accounts for the important 
physics. Furthermore, the various phase shifts between signals from 
different probes match in the one- and two-cell cases. 

5.6 Comparison with experiments 


Experiments [Tadjfar & Others, 1985] at Re k = 310 and k/6* = 
1.6, with an array of spherical roughness elements similar to the 
standard density case used in this study, reveal Stk = fk/Uk = 0.08 
using the undisturbed Blasius velocity and 0.23 using the measured 
velocity at the top of a sphere. These results compare favorably 
with Stk = 0.23 with isolated spheres [Mochizuki, 1961], but dif- 
fer from the Stk = 0.33 numerical results. This disparity cannot 
be resolved without additional experimental results using a different 
roughness density. The numerical and experimental results agree in 
one area, however. The peak amplitude occurred near 1.5 k from the 
wall in both cases, although Tadjfar took no measurements within 
the roughness elements. 

6 Conclusions 

Since identical oscillation frequencies were obtained with peri- 
odic boundary conditions applied across both one- and two-cell do- 
mains, the domain length does not govern the oscillation frequency; 
instead, the primary instability associated with a periodic array of 
spheres protruding from a flat surface occurs at a frequency which 
depends on the streamwise spacing between elements and the lo- 
cal velocity at the elevation of the sphere’s crest. This also means 
that the characteristic frequency remains independent of the isolated 
sphere shedding frequency, and any similarities between the two are 
merely coincidental. 
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This oscillation frequency, called the passing frequency for the 
time required for a fluid particle to pass from one sphere to the next, 
may also be thought of as a forcing frequency due to disturbances 
generated in the wake region behind the individual spheres. The 
response where the oscillation frequency equals the forcing frequency 
is called a lock-in. In the lock-in mode, the frequency of oscillation 
must lie near the shedding frequency of a similar body located in an 
infinite domain — in this case, the similar body equals a single sphere 
attached to a flat surface. 

At Rek = 175, no streamwise vortices based on the flow’s inertia 
exist; instead, the flow field contains regions of streamwise vorticity 
explained entirely by continuity considerations. 
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